FIG. 1: Tiling of the square lattice with 13-site clusters. The 
NN hoppings within a cluster are represented by full lines, 
and inter-cluster hoppings by dashed lines. Of course, much 
simpler tilings (e.g. square) are possible. 

cluster is formally the quotient j/T, We will use the nota- 
tion r and R for the lattice sites of 7 and T, respectively, 
and L will denote the number of lattice sites within a 
cluster. An example of clustering of the two-dimensional 
square lattice with L = 13 is shown on Fig. @. Tiling the 
lattice with rectangular clusters is obviously simpler, but 
being able to compare the results of different tilings is 
important, in order to separate real physical effects from 
artefacts caused by a particular cluster shape. Each lat- 
tice site r of a cluster can also host a number r£ orD > 1 
of orbitals (as in a many-band model) , although the case 
"■orb = 1 will be most frequent. 

The class of Hamiltonians that can be treated by this 
approach has the form H = Hq + V , where 

h» = J2 h * . v = E ^ " ' n. ( " '' (i) 

R R,R' 

a . b 

Here a, b label different electron orbitals (or sites) within 
a cluster, and range from 1 to n orD £. is a Hub- 
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bard Hamiltonian defined on a single cluster and V is a 
hopping term between clusters, V^jj 11 being the hopping 
amplitude between orbital a of cluster R and orbital b of 
cluster R'. Hopping can be of any range, and multi-band 
models are treated on the same footing as the one-band 
model: The indices a, b may refer to different lattice sites 
on a cluster, or to different bands (orbitals) within a clus- 
ter, or both. In order to stay as general as possible, we 
will call them "orbital indices" . Spin indices are implicit 
in hopping terms and suppressed where appropriate in 
order to lighten the notation. 

In the simplest case, treated explicitly in Refs |]]^, the 
cluster consists of a single site (7 = T and n or b = 1) and 
the perturbation V is a nearest-neighbor hopping : 



yr 



-t (r,r'N.N.), (2) 



In that case, the perturbative treatment of V coincides 
with the so-called strong-coupling perturbation theory 
(SCPT) of the Hubbard model and has been extensively 
studied in Ref. || Additional references on the strong- 
coupling perturbation theory, in particular earlier at- 
tempts at its systematization, can be found therein. We 
review the main results in the following subsection. In 
practice, SCPT provides approximate expressions for the 
various electron Green functions in terms of the hop- 
ping integral and of various atomic Green functions (i.e., 
Green functions of the one-site Hubbard model). The 



general results of SCPT can be immediately applied to 
the class ([]]) of Hamiltonians, with the difference that 
atomic Green functions are replaced with exact Green 
functions of a Hubbard model defined on a small cluster. 



In Sect. Ill, we make a series of important remarks on 



the method. In Sect. IV we explain how the Green func- 
tion on the cluster is calculated numerically. We com- 
ment on multi-band models in Sect. We conclude in 



Sect. VI with a calculation of integrated quantities, like 
the energy density and double occupancy. 



II. DERIVATION OF THE METHOD 
A. The strong-coupling expansion 

In this section we review the main results of Ref. 
on strong-coupling perturbation theory, in a notation 
adapted to the present situation. It is assumed that we 
are dealing with a Hamiltonian of the form (Q) . Our goal 
is to place the approximation made in CPT in the con- 
text of a systematic strong-coupling perturbation theory, 
i.e. to explain in what way it constitutes a leading order 
approximation. 

In the path-integral formalism, the partition function 
is expressed as a functional integral over Grassmann vari- 
ables 7R, a that correspond to the creation operator cr : 



Z = f [d 7 *d 7 ] exp - J dr J 7R» (d T - n) 7Ra(r) + £ H R ( lRallRa ) + £ V**' lRalR , b \ 

J J ° ^R,a R R,R' j 



(3) 



In order to lighten the notation, we use greek indices to denote sets such as (R, a, r), and use bra-ket notation. For 
instance, 



I *r E KT'7r>)7R' 6 (t) =J2 V7> = (7M7) 

•'O D D' 11.1/ 



(4) 



R.R 

a , 6 



Strong-coupling perturbation theory, as formulated in Ref. |^J^, proceeds first by a so-called Grassmannian Hubbard- 
Stratonovich transformation, which amounts to expressing the perturbation (7IVI7) as the result of a Gaussian integral 
over auxiliary Grassmann fields ip-R a (t) , ip Ra (t) : 



e -<7|v| 7 > = det v j [ dr p* d ^] exp [{iplV-^ip) + (V>|7) + (7IV)] (5) 

In terms of the auxiliary field, the partition function becomes, up to a normalization factor, 

Z = Z Q J [d?P*diP}e wv ~ lw ^W+frl*^ (6) 

where (• • -)o means an unperturbed average over the original electron field 7. Denoting- by (• • -)o, c the cumulant 
averages, and owing to the block-diagonality of H°, the above average can be rewritten astl: 

00 ^ _g R 

fl=l R{a l: a[} J0 1=1 



FIG. 2: Diagrams associated with the terms of order V° , 
V 1 (above) and V 2 (below) in the self-energy 3 of auxiliary 
fermions, in strong-coupling perturbation theory. Here only 
the interaction terms R — 1,2, 3 of Eq. (0) contribute. The 
term R — 4 (i.e., G?' 4 ') would come in at order V 4 . 

Translation invariance along the superlattice T allows 
us to express V and Q in terms of a wavevector Q be- 
longing to the reduced Brillouin zone (BZr, the Brillouin 
zone of the superlattice) instead of cluster indices. In 
particular, the hopping term may be expressed as 

^(Q)=EC eiQ,R ( 12 ) 

R 

where a superlattice site '0' has been chosen as the origin. 
In this representation, Eq. (|l0|) becomes 

g a ,m^)= ( — (is) 

This is the starting formula of Cluster Perturbation The- 
ory on homogeneous systems. 

B. Lattices, superlattices and wavevectors 

The Green function Ga,b(Q,z) is in a mixed repre- 
sentation: direct space within a cluster and reciprocal 
(Fourier) space between clusters. A pure Fourier repre- 
sentation is preferable, in terms of wavevectors belonging 
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to the Brillouin zone BZ 7 of the original lattice. However, 
the perturbative treatment breaks translation invariance 
on 7, because it singles out hopping terms between clus- 
ters, while translation invariance over the superlattice T 
is preserved. As a result, the Green function Ga,b(Q>z) 
will depend on two wavevectors k and k' of the Brillouin 
zone, except that k — k' must belong to the reciprocal 
superlattice T* . We will demonstrate the following (the 
frequency argument is suppressed here): 



S(k,k') 



1 



L L 
s=l a,b=l 



6(k-k' + d s )g ab {k)e 



-ik r a ik' r b 



(14) 

where (i) q s belongs both to the reciprocal superlattice 
r* and to the original Brillouin zone BZ 7 (L possible 
values) and (ii) r is the position of the lattice site as- 
sociated with the orbital index a. For simplicity, we will 
assume a one-band model, so that the orbital index a 
is also a site index. This assumption will be relaxed in 
Sect. below. 

The demonstration of Eq. (|l4]) is straightforward. We 
simply apply the basic Fourier transform 

1 



CRa = 



,c(kK 



k-(R+r a ) 



(15) 



where N is the number of sites on the superlattice (iV 
oo). Applied to the Green function, this yields 

NT ^ 



S(k,k') 
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R,R' 

a , b 



/->R.R 
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ik'-(R'+r b ) 
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<S(K - Q)<S(K' - Q) 



a.b 



5ab(K) C 
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ik.'. 



^(K-K') 



(16) 



where the wavevector k has been decomposed in a unique 
fashion as k = K + k, where K belongs to the re- 
duced Brillouin zone BZp and k belongs to the recip- 
rocal superlattice T*, and likewise for k'. Note however 
that V(K.) = V(k), since the definition ( |l2| ) is invari- 
ant with respect to translations of Q by a an element 
of the reciprocal superlattice, such as k — K; therefore, 
Gab(K) = ^ab(k). Since 



5(K-K')=^S(k-k' + q s ) 



(17) 



we recover Formula (|l4|). 

The q s = term in the above sum is the CPT approx- 
imation to the translation-invariant Green function: 



1 L 

S C P T (k,z) = 7 E 5ab(Kz)e-*< r «-^ 



(18) 



a,fc=l 



This, together with Eq. jl^), is the central formula of 
Cluster Perturbation Theory. In the remainder of this 
paper, we will drop the index 'CPT' on the Green func- 
tion or related quantities. 

In pratice, we are more often interested in the spectral 
function 

A(k,w) = -2 lim Im0(k,u; + M7 + /z) (19) 

where \i is the chemical potential. Unless particle-hole 
symmetry is present, the chemical potential is not known, 
since the calculation is done at fixed filling. It must be 
calculated from the density of states (this will be ex- 
plained in Sect. VI below). For this reason, it is more 
pratical to deal with the shifted spectral function 

A(k, u)) = A(k, u> - £i) = -2 lim lmg(k,oj + irj) (20) 

which does not involve /x in any way (Note that in the 
present paper, the Hamiltonian excludes the chemical po- 
tential term as well). 

The q s 7^ terms are not considered here, but they 
are a measure of the breaking of translation invariance 
caused by the different treatment given to intercluster 
and intracluster hopping. It has been verified that the 
spectral weight associated with those terms is nonposi- 
tive, and that it integrates to zero. 



III. GENERAL REMARKS 

A. Limiting cases 

The CPT Green function ( p"8| ) is exact both in the 
strong and weak-coupling limits. This may look para- 
doxical at first, since it is a perturbative result that is 
expected to be exact only when the hopping \ ' goes to 
zero. However, if we set U = from the start and use 
ordinary 'weak-coupling', perturbation theory (i.e. treat 
inter-cluster hopping as a perturbation of in-cluster hop- 
ping), then obviously V is the exact electron self-energy 
and the basic relation Q^ 1 = G^ 1 — V follows. In other 
words, and using the notation of a nearest-neighbor Hub- 
bard model, the previous relation may be viewed either 
as resulting from the lowest order term (5 = G) in a t/U 
expansion of the auxiliary fermion self-energy or from 
the exact U = electron self energy (E = V). Need- 
less to say, verifying that the method gives the known 
U = spectrum is a practical test that is often carried 
in practise. Incidently, Dynamical Mean Field Theory 
(DMFT)Q also is exact in those two limits. The relation 
between the two approaches (DMFT and CPT) remains 
an open question. 



B. Nature of the approximation 

It is somewhat difficult to qualify the degree of approx- 
imation achieved by Cluster Perturbation Theory. The 



formula (|1C|). This is complicated by the presence of the 
same physical parameter (the hopping t) both within the 
cluster (in which it is treated exactly) and between the 
cluster (where it is treated perturbatively) . In principle, 
the CPT could be systematically improved by taking into 
account the higher-order diagrams of Fig. § but we have 
already pointed out that their exact numerical calcula- 
tion is impractical. In practice, the CPT approximation 
is controlled most effectively by the cluster size L. 

A look at Fig. | is enough to demonstrate the tremen- 
dous advantage of treating clusters exactly. Here we 
treated the one-dimensional, nearest-neighbor Hubbard 
model at half-filling, with clusters of sizes 1 and 2. If 
L = 1, one falls back to the zeroth-order approximation 
of strong-coupling perturbation theory; this leading order 
coincides with the Hubbard-I approximation!] : 



uj - 2t cos k - U 2 /Auj 



(21) 



FIG. 4: Spectral function of the two-dimensional Hubbard 
model (U = 8t, t' = -OAt, n = 2/3) as calculated by CPT on 
four different clusters. The momentum scans are indicated. 



The difference between this approximation and the 
atomic limit lies in the dispersion acquired by the two 
Hubbard bands, and in the momentum-dependent trans- 
fer of spectral weight between the two bands. Note that 
this is also the zeroth-order case of a calculation that was 
pushed to fifth order (i.e. {t/Uf) in Ref. |. If we now 
consider the spectral weight obtained from a cluster of 
size L = 2, we already see the main features of the pre- 
sumed exact result, in particular the weak shadow bands 
that disperse completely from to 7r, i.e., on the reduced 
(magnetic) Brillouin zone. Of course, the superexchange 
J = At 2 /U emerges from the exact solution of the two- 
site cluster and has an impact on the corresponding CPT 
spectrum. The L = 2 spectrum is also closer to the pre- 
sumed exact result than that obtained by strong- coupling 
perturbation theory at order (t/U) 3 (Fig. 3 of Ref. l|). 

In Fig. ||, we show the spectral function of the two- 
dimensional Hubbard model (with NN and NNN hop- 
ping and density n = 2/3) for four different clusters (il- 
lustrated on the upper right corner of each plot). Minor 
details may vary from cluster to cluster, but the gross 
features stay the same : (i) the general dispersion of the 
main band; (ii) the two closely separated bands between 
the points T and X, the upper one very close to the 
Fermi level; (iii) the two well-separated (Aoj ~ U) fea- 
tures at the antifoerromagnetic wavevector M = (tt, 7t). 
The larger the cluster, the more poles contribute to the 
spectral function, making it smoother. 



C. The question of boundary conditions 



approach is exact in the limits U/t — 0, t/U — and 
L — > oo. The perturbation parameter is th e interclus- 
ter hopping V and, according to Sect. [I A above, the 



V dependence has been dropped from the auxiliary field 
self-energy S, since only the first diagram of Fig. || has 
been kept. However, the CPT Green function itself con- 
tains terms of all orders in V, as seen from the basic 



In its usual formulation, Cluster Perturbation Theory 
requires that the cluster Green function be calculated 
with open boundary conditions, as opposed to periodic 
boundary conditions. However, this makes the numer- 
ics a little more demanding, since the wavevector within 
a cluster is no longer a conserved quantum number and 
the Hilbert space to work with is larger by a factor L. 



FIG. 6: Density plots of the spectral function of the two- 
dimensional Hubbard model with NNN hopping t' — — O.At, in 
the first quadrant of the Brillouin zone. In (a), the density is 
n = 8/9 (i.e., 11% hole doping) and U = lit. In (b), n = 10/9 
(i.e., 11% electron doping) and U = At. In both cases, the 
noninteracting Fermi surface is shown (dashed curve). The 
spectral function was integrated over frequency in a interval 
A = below the Fermi level. 



hopping t' — — OAt and n = 8/9. This figure should he 
compared with ARPES results on Ca2-zNa ;r CuC)2Cl2Ll, 
a hole-doped material. One notices a disappearance of 
the Fermi surface in the antinodal directions (tt, 0) and 
(0, 7r), which coincides with the opening of a pseudogap 
in those directions. By contrast, Fig. shows the same 
spectral function for n — 10/9 and U = At. There, the 
Fermi surface disappears first on "hot spots" lying on 
the (tt, 0) — (0, 7r) line; at higher values of U, the Fermi 
surface disappears in the nodal direction, but remains 
in the antinodal direction, in contrast to the hole-doped 
case. Tbis should be compared with ARPES results on 

nccoK 



E. The t — J model 

The derivation of CPT given above is specific to 
Hubbard- like models, i.e., models with only on-site re- 
pulsion and hopping terms. The t — J model, because 
it involves correlated hopping, seems excluded from the 
method. However, CPT has been applied with some suc- 
cess to the t — J model in Refs (y],[L2[ and this raises the 
question as to why it can be successful in this case. No 
rigorous answer can be given. However, one may first 
remark that, near the line J = At 2 /U, the t — J model is 
a low-energy limit of the Hubbard model, and therefore 
CPT wouls be expected to work, at least at low energy. 
Second, in CPT, the correlated hopping of the t— J model 
is replaced by ordinary hopping only between clusters, 
and for a large cluster this may be of little consequence. 
Finally, it may very well be that the basic formula ( |l3| ) 
is also the lowest-order approximant to a systematic ap- 
proximation scheme to the t — J model unknown so far, 
albeit different from the Hubbard-model strong-coupling 
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perturbation theory. 



F. Finite temperature 

In principle, CPT can be used at finite temperature. 
Ground state averages are simply replaced by thermal av- 
erages (in the canonical or grand canonical ensembles). 
If the cluster is small enough that all eigenstates can 
be calculated, then the finite-temperature extension is 
a straightforward matter. In those cases, computation 
time is already short, and a finite temperature is not too 
much of a numerical burden. For instance, the entropy 
of the two-dimensional Hubbard model and related ther- 
modynamic quantities can be reasonably well calculated 
on a 2 x 2 clusterti If the cluster is large enough that 
the Lanczos algorithm is necessary, then the same tech- 
niques used in exact diagonalizations at finite tempera- 
ture must be usedtj. These imply a sampling of initial 
Lanczos states. In that case computation time may be- 
come very important, and some form of parallelization 
necessary. Work on this topic is in progress. 

Incidently, CPT is independent of the way the cluster 
Green function is calculated. Exact diagonalizations us- 
ing the Lanczos algorithm are appropriate at zero tem- 
perature, but another technique, like Quantum Monte 
Carlo (QMC), could be used at finite temperature. How- 
ever, the Green function obtained from QMC is defined in 
imaginary time, and the continuation to real frequency 
(via the maximum entropy method) is not only time- 
consuming, but has problems resolving finer structures 
in the spectral function. At this point it is not clear 
whether a Monte-Carlo driven CPT would be a produc- 
tive method. 



IV. THE CLUSTER GREEN FUNCTION 

A. The Lanczos method 

In this subsection, we will briefly review how to cal- 
culate the cluster Green function G a b(uj) numerically. 
First, since the Hamiltonian conserves separately jVj and 
iVj_, the numbers of spin- up and spin-down electrons on 
the cluster, there is a practical advantage in calculating 
the Green function at fixed values of and iVj, : it 
saves memory (these numbers are not longer conserved 
once perturbation theory is applied, because of the inter- 
cluster hopping term) . One must first identify the sector 
that contains the ground state. In the absence of mag- 
netic field, the ground state is expected to be a singlet, 
i.e., to lie in the = N± sector. An interesting-ex- 
ception to this rule is indicated by Lieb's theoremlij : 
on a bipartite lattice with A sites on the first sublat- 
tice and B sites on the second sublattice, the ground 
state of the half-filled, repulsive, nearest-neighbor Hub- 
bard model has spin \A — B\, 



The second step is to use the Lanczos algorithm (or the 
equivalent) to calculate the ground state in the cho- 
sen sector. One needs to implement a practical encoding 
of the states. The Hilbert space can be factorized into 
spin-up and spin-down parts : V = <S> Vj,. A cluster of 
L sites with N- 



d = 



N has dimension 



N\{L-N)\ 



(22) 



For instance, a half-filled 12-site cluster has dimension 
d = 853 776, and storing one state in memory (in dou- 
ble precision) requires 6.5 MB. At 16 sites, this figure 
becomes 1.2 GB, which is impractical unless additional 
symmetries are used. However, since the method usually 
demands open boundary conditions, translational invari- 
ance is of no use, and rotational or inversion symmetry 
are typically too much trouble to implement for what 
they save in memory. Once a good approximation to |fi) 
has been obtained, one may then proceed to calculate the 
electron and hole parts of the Green functions : 



Gab(z) 

G h ab (z) 



= G e ab (z) + G h ab (z) 



= (n\4 



z-H + E 
1 

z + H-E a 



c a \n) 



(23) 



where Eq is the ground state energy and z a complex 
frequency. In practice, this is done as follows : one con- 
structs the states \<j>) = c b \fl) and \<j>') — <4|f2). Then 
\4>), after being normalized, is used as an initial state in 
a second Lanczos iteration, which produces a sequence 
of orthonormal states \(j> m ), in terms of which the Hamil- 
tonian H is a tridiagonal matrix. At each Lanczos step, 
the product (4>'\4> m ) is calculated, and the states \<j) m ) are 
not stored, except for the two most recent ones, which are 
used to calculate the next state. After a sufficient num- 
ber of iterations (typically a few hundred) , one may write 
with good accuracy 



G e ab {z) 



1 



'Z-H + En 



(24) 

Since |^>o) oc \<fi), one then has, for a given value of z, 
to solve a tridiagonal system of equations to find the 
column- vector {x m }, and this is relatively efficient nu- 
merically. The same procedure is repeated for the hole 
part G a ] b (z). In practice, one stores in memory the tridi- 
agonal matrix ((j) m \H\(f) n ) for each pair of sites (a, 6), as 
well as the products (<j>'\(j) m ). This means that the Lanc- 
zos procedure is used only at the beginning, before the 
sweep over frequencies is performed. 

Note that the full Hamiltonian can easily be stored in 
memory during the Lanczos procedure. Indeed, it can be 
expressed as 



H = K\ ® 1 + 1 ® K i + Vboui. 



(25) 



FIG. 7: Chemical potential as a function of density in the 
one- dimensional Hubbard model, as calculated by CPT using 
a supercluster of ten ten-site clusters (dots). This is to be 
compared to the exact result obtained by differentiating the 
exact ground-state energy density with respect to density, fol- 
lowing the prescription of Ref. U & (solid line). 



reads 

-, n orb L 

£$ T (k, z) = - £ £f„(k, z)e-^ r — ) (26) 

a,b=l 

from which we derive the many-band spectral function 
A afl (k,uj) = -2 lim Im£ Q/3 (k,Lj + ^ + a) (27) 

If an electron is annihilated in a combination of bands by 
the operator u a c a (k), then the relevant spectral function 
is (summation over repeated indices is implicit) 

p a /3A af3 (k,U)) paf3 = U a U*p (28) 

More generally, p a p can be any suitable density matrix 
(with unit trace) and may correspond to a pure state 
(like above) or to a mixed state, depending on the phys- 
ical process under study. Giving it a suitable momen- 
tum dependence may also be a way of including matrix- 
element effects, as in Ref. But the simplest case would 
be that of a photoemission process in which the photo- 
electron comes from any band with equal probability, in 
which case the density matrix is unity. This scheme was 
applied in Ref. ^ (albeit using periodic boundary condi- 
tions on the cluster) in studying the spectral function of 
Sr 2 Cu0 2 Cl 2 . 

We illustrate in Fig. || the spectral function for a three- 
band Hubbard model. The band parameters are those 



proposed in Ref. 19: (i) A hopping £ = 1.6 eV between 
Cu and O orbitals; (ii) a next-nearest-neighbor hopping 
t' = 1.1 cV between O orbitals, both diagonally and 
across a Cu atom; (hi) a shift e p — —3 eV of the O 
bands with respect to the Cu band. In addition, we put 
an on-site repulsion U = 4 eV on the Cu orbital only. 
This is meant to be a realistic model of Cu02 planes in 



FIG. 8: Spectral function for a three-band model of CuC>2 
planes. The wavevector scan is the same as usual: V = (0, 0), 
X — (tt,0) and M — (tt,tv). Exceptionally, we plotted the 
shifted spectral function A(k, uS) = A(k, U) — fj.) (the value of 
/x is indicated). The U = dispersion relations are indicated 
(thick lines). 



YBCO, except for the value of U, which we set arbitrar- 
ily. Filling was set to 11/12 (i.e. n — 11/6) for the three 
bands, which means, in a free electron language, that 
the upper band has filling 3/4 and the lower two bands 
are completely filled. We have also plotted the noninter- 
acting (U — 0) dispersion relations for the three bands 
(thick lines). In the absence of interaction, the spectrum 
consists, for each wavevector, of three equal-amplitude 
delta- function peaks, one for each band. One sees that 
the interaction has affected mostly the upper band, but 
that the lower (mostly oxygen) bands are also affected. 
In particular, the upper band shows extended spectral 
weight from T to M, and a gap of order ~ U/2 to a 
narrow, almost dispersionless feature between M and X, 
lying just above the Fermi level. The lowest band shows 
a small gap near its inflection point (M), and some por- 
tions of the lower bands are somewhat raised by the in- 
teraction. This calculation is based on a four-site cluster 
with three orbitals per site, 25 hopping terms within the 
cluster and 30 with adjacent clusters. 



VI. ENERGY, DOUBLE OCCUPANCY AND 
CHEMICAL POTENTIAL 



In this section we will explain how to calculate the 
average energy per site and double occupancy from the 
one-particle Green function, at zero and finite tempera- 
ture. The basic ingredient of this calculation is the rela- 
tion between the ground state eusrgy density e and the 
zero-temperature Green functioned: 



dk duj 

{2n) d 2tt 



e*"° (w + £ k )G T (w + iO + sign(w-/i)) 



(29) 



(we assume here that G-\ — G\). Since we use a 
complex- frequency Green function throughout, the usual 
(time-ordered) real- frequency Green function is + 
«0 + sign(tj — ju)), whereas the retarded Green function is 
simply G" T ct (w) = G T (u + i0+). With the help of the 
spectral representation 



G(k,z) 



du> A(k,uj) 



2tt z — u> — /i 

we find, after integrating around the upper half-plane, 
dk f° duj 



(30) 



(2ir) d 



2vr 



(lu + fi + e k )A(k,Lu) 



(31) 



The kinetic energy per site £kin may be found in a similar 
way: starting from 

f dk 

£ki„ = 2 J _- ek (Q|4 iT c kiT |n) (32) 

(the factor of 2 comes from the sum over spin), we use 
the definition of the time-ordered Green function 



G a (k,t) = -i(fi|Tck, CT (0)c kiCT (t)|O) 
to express it as 



(33) 



£kin = 2i 



dk 

(2^ 
dk 

W) 

Therefore, by subtracting ([m]) from ([n]), one finds the 
potential energy per site in terms of the spectral function: 



2tt 1 v 



«0 + sign(cj — fj,)) 
(34) 



-pot 



dk 

(2tt) c 



2^ 



(u + [i - e k )A(k, u) (35) 



and the double occupancy is simply V = e po t/U. 

Before calculating the energy and the double occu- 
pancy from formulas ( |3l"| , |3"5| ) , the chemical potential fi 
must be known. But except for special cases with 
particle-hole symmetry (half-filling, with no NNN hop- 
ping), \i must be calculated from the electron density n 
by imposing the constraint 



n = 2 



dk 

{2n) d 



(36) 



where A(k, u>) is the /i-independent, shifted spectral func- 
tion of Eq. (p0[). In practice, this is done by numerical 
integration, where the momentum integration is carried 
before the frequency integration, since the evaluation of 
the Green function is much faster if momentum sweeps 
are performed at fixed frequency. 

The formulas ( |3l| , p5| ) can then be applied to the CPT 
spectral function A(k,uj). One basically needs to com- 
pute the two integrals 



dk 

(2-K) d 

dk 
(27r) d 



2^ 



e k A(k,iv) 



did 

—wA{k,uj) 

Z7T 



(37) 



FIG. 10: Same as Fi g. H . this time for the double occupancy, 
calculated from Eq. (p5^. An extrapolation of the results to 
infinite cluster size (L — > oo) using a quadratic fit in terms 
of 1/L is also shown, and is accurate to within 0.5%. exact 
diagonalization results for L = 8 and L = 12 are also shown. 



All these integrals depend on the line broadening r), 
and are increasingly difficult to compute as -q decreases. 
Extrapolating towards 77 = from moderate values of 
7] (between 0.1 to 0.02) is the method of choice. These 
tricks are also used in the prior calculation of the chemical 
potential from Eq. ([56]). 

Fig. ^| shows the relative difference between the com- 
puted ground-state energy density of the one-dimensional 
Hubbard model at half-filling and the known exact result 
computed from Ref. [f(| Even with as few as 4 sites per 
cluster, an accuracy of 1% is achieved. Clusters of 12 
sites bring this down to 0.2% or less. Notice that the 
energy density of a bare cluster (no intercluster hopping) 
with periodic boundary conditions (to eliminate edge ef- 
fects) is not nearly as close (dashed lines) to the exact 
result as the CPT prediction. Fig. [l(] shows the same 
calculation for the double occupancy. Here the accuracy 
obtained is less impressive (4% for a 12-site cluster), but 
a simple quadratic extrapolation as a function of inverse 
cluster size yields the exact result to within 0.5%. More- 
over, the CPT predictions converge to the exact result in 
a smoother fashion (as a function of t) than exact diago- 
nalization data. Cluster Perturbation Theory is therefore 
a quite accurate way of computing global quantities like 
the energy density or the double occupancy. 



VII. CONCLUSION 

Cluster Perturbation Theory is an economical method 
for calculating the approximate spectral function of Hub- 
bard models. Its main advantage over exact diagonaliza- 
tions is the access to a continuum of wavevectors and a 
better approximation to the thermodynamic limit. This, 
for instance, allows for an investigation of pseudogap phe- 
nomena (Fig. ^). It is also much faster than Quantum 
Monte Carlo calculations at very low temperatures and 
may be formulated directly in terms of real frequencies. 
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Its accuracy has been demonstrated here in the calcula- 
tion of the ground-state energy of the one-dimensional 
Hubbard model. In our opinion, it is a method of choice 
for exploring parameter space and trying to connect real 
materials to Hamiltonians with local interactions. 
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